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SUMMARY 

The flux-integral method is a procedure for constructing an explicit, single-step, forward-in- 
time, conservative, control-volume update of the unsteady, multidimensional convection- 
diffusion equation. The convective-plus-diffusive flux at each face of a control- volume cell 
is estimated by integrating the transported variable and its face-normal derivative over the 
volume swept out by the convecting velocity field. This yields a unique description of the 
fluxes, whereas other conservative methods rely on nonunique, arbitrary pseudoflux- 
difference splitting procedures. The accuracy of the resulting scheme depends on the form 
of the sub-cell interpolation assumed, given cell-average data. Cellwise constant behaviour 
results in a (very artificially diffusive) first-order convection scheme. Second-order 
convection-diffusion schemes correspond to cellwise linear (or bilinear) sub-cell 
interpolation. Cellwise quadratic sub-cell interpolants generate a highly accurate convection- 
diffusion scheme with excellent phase accuracy. Under constant-coefficient conditions, this 
is a uniformly third-order polynomial interpolation algorithm (UTOPIA). 


THE FLUX INTEGRAL 

Consider the cell-average value of the transported scalar, <t > , at a reference (central) cell, C. 
In two dimensions, an exact , single-step, explicit update can be written for the “new” 
(superscript +) cell value: 



<t>c + FLUX w (zj) - FLUX w (/+ l,y) + FLUX ^, 7 ) - FLUX f (iJ+l) C 1 ) 


using standard index and compass-point notation. Note that this is strictly conservative in 
that the east-face convective-plus-diffusive flux of cell C, at (zj), is identical to the west- 
face flux at ( 1 + 1 , 7 ); similarly for the north- and south-face fluxes. In (1), the west-face 
flux, for example, is given by 

FLUX.OV) = (c„*.> - h («. (||) w ) (2) 


where the angle-brackets represent time-averages over At, and, assuming (for convenience) 
a uniform square mesh of side h, the west-face normal-component Courant number is 

c = (3) 


and the west-face nondimensional diffusion parameter is written in terms of the (scalar) 
diffusivity, D w , as 


a 


w 


DJt) At 

h 2 


(4) 


with analogous definitions for the south face. 

The convective contribution in (2) is equivalent to the total “mass” of <f> swept 
through the west face along particle paths (or streamlines, in steady flow) over At. In 
principle, one could trace the particle paths backwards to the earlier time-level, for each 
face. This is shown, schematically, in Figure 1. Then the (exact) purely convective 
contribution is equivalent to integrating 4>(x,y) at the earlier time-level over the area (or 
volume, in three dimensions) swept out by the particle paths: 

= J J 4>(x,y) d (x/h) d (y/h) (5) 

PPA 


where PPA stands for particle-path area. 
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The flux-integral method now approximates (5) by replacing the particle-path area by 
the flux -integral parallelogram (FIP) by assuming the convecting velocity field to be locally 
constant (in both space and time) in the vicinity of the face in question. This is shown in 
Figure 2; note that the parallelogram is defined by the local (space-time-averaged) Courant 
number components, c„ and c^ (taken as both positive in the case shown). The flux-integral 
convective approximation is thus 

(c^ <t> w ) * J J 4>(x,y) d (x/h) d(y/h) (6) 

FIP 

A similar approximation for the diffusive contribution results in 

- * («- QgL) " - % J l % « (7) 

w FIP 

where a w is an appropriate average. If the sub-cell behaviour at the earlier time-level, 
4>(x,y), were known in complete detail, (6) would represent an exact flux due to pure 
convection by a constant velocity field. For non-zero diffusion, it turns out that (7) 
represents the diffusive flux to third order— provided the sub-cell behaviour is known in 
enough detail. The major task in the flux-integral method is thus an interpolation problem: 


Given a set of cell-average values, 
estimate sub-cell behaviour in an accurate 
(and, ideally, shape-preserving) manner, 
while observing the cell-average constraint: 


4>(x,y) d (x/h) d (y/h) 

cell 


<f> for all cells 


with an analogous formula in three dimensions. 


( 8 ) 


For constant-density flow, the face-normal component Courant numbers used in 
constructing the flux-integral parallelograms should satisfy a discrete continuity equation for 
each cell: 

cjtj) - Cji+Uj) * cjfj) - cjij+l) = 0 (9) 
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Since the area of a flux-integral parallelogram is proportional to the face-normal Courant 
number component (e.g., the shaded area in Figure 2 is c m h 2 ), an initially constant scalar 
field will remain constant everywhere (to machine accuracy). This can be seen in Figure 3, 
where the sum of the “inflow areas” equals the sum of the “outflow areas” (irrespective 
of the local individual face-transverse Courant number components). 

In the following sections, a number of different sub-cell interpolants are explored. 
A cellwise constant interpolant (locally equal to the cell average) results in a (very 
artificially diffusive) first-order convection scheme; modelled physical diffusion is absent, 
to a consistent order. This is not a viable scheme for practical CFD calculations. But, 
because of its simplicity, it is instructive to explore the flux-integral method in this case. 
Bilinear downwind - weighted sub-cell interpolation results in a two-dimensional analogue of 
the Lax-Wendroff scheme [1]. A cellwise quadratic interpolant over each cell generates a 
convection-diffusion scheme that is formally third-order accurate under constant-coefficient 
conditions. These three methods are briefly compared using the well-known “rotating hill” 
test problem. 


FIRST-ORDER CONVECTION 

Fluxes will be calculated for the west face of cell C. Entirely analogous fluxes for the south 
face can be written down using appropriate (x,y) permutations. Unless otherwise noted, the 
Courant number components will be taken as both positive. Referring to Figure 4, the west- 
face convective flux integral corresponding to (6) is seen to consist of three parts 

FLUX w (t j) = - I 2 * I 3 (10) 

where /, is the integral over the rectangular area, 1246, in cell W; 1 2 is over the triangular 
area, 123, in cell W; and I 3 is over a similar triangular area, 456, in cell SW. Assuming <j> 
to be cellwise constant (equal to the respective cell average), the integrals in (1) are 
proportional to the respective areas times the local cell-average value. This gives 
FIRST-ORDER FIM: 

FLUX„(/j) = e_ ^ < u > 

or, on rearrangement, 
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FLU X„(iJ) = c„[k- c -f (K - *„)] 


( 12 ) 


The term in square brackets can be considered to be the effective average convected face 
value. Note that it consists of the one-dimensional, first-order upwind (“donor-cell”) value, 
4> w , modified by “a transverse-gradient term proportional to the transverse Courant number 
component at the face. It should be clear how the formula changes for other combinations 
of signs of c m and c^. To a consistent order, there is no (physical) diffusive flux, since 
cellwise constant behaviour implies a zero-gradient within each cell. 

Substituting (12), and the analogous formula for FLUX,, into (1) gives an overall 
(constant-coefficient) convective update equation: 


This is identical to a semi-Lagrangian update [2], using bilinear interpolation around the 
departure point, collocated at node values: <f> c , 4> w , ^ , and <f> s (located at the centroids 
of the respective cells). [For cellwise constant, linear, or bilinear interpolants, node values 
are equal to the respective cell-average values. This is not the case for higher-order 
interpolants.] 

Using an appropriate upwinding strategy for other convecting velocity directions, it 
is not hard to show that the von Neumann stability condition for this scheme is given by a 
square region in the (c x , c y ) plane: 


SECOND-ORDER METHODS 

Second-order convection-diffusion methods result from assuming cellwise linear or bilinear 
sub-cell behaviour. In this case, it is convenient to introduce local, normalised coordinates 


4>C = 4>C " C x (&C ~ 4>w) ~ C y (^c ~ ^s) 

+ C x C y (h ~t W -4> S + 4 >Sw) 


(13) 


|cj <1 and | C y | < 1 


(14) 


in cell W: 


f i - (■• - 1) 


(15) 


and 



( 16 ) 
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as shown in Figure 5. Note that the central cell, C, is located at (i,j). The top of the flux- 
integral parallelogram is represented by 

*!„«) = 0.5 ♦ fs (f - 0.5) (17) 

^IW 

A general bilinear sub-cell interpolant within cell W takes the form 

HS,n) = C x + C 2 £ + C 3 7, + C 4 £77 (18) 

The cell-average constraint, (8), implies 

c, - K (19) 

The slope-constants, C 2 and C 3 , and the twist-constant, C 4 , can be chosen in a number of 
different ways. For example, a two-dimensional analogue of Fromm’s method [3] results 
from choosing 

c 2 - i (*c - o < 20 > 



and 

C 4 = 0 (22) 

Note that this represents a symmetrical distribution with respect to cell W, independent of 
velocity direction. Upwind or downwind weighting can also be used. In the interest of 
brevity, only one scheme will be considered here in detail. This is based on downwind- 
weighted bilinear interpolation. For example, throughout cell W, the interpolant is 
collocated (for c w , > 0) at node values: <f> c , 4> N , and (in addition to <f> w ). This 
turns out to be a two-dimensional generalisation of the Lax-Wendroff method. For positive 
Courant number components, the interpolant within cell W is 
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( 23 ) 


<t>(£,v) = 4> w + ( <f> c ~ 4> w ) £ + - <l> w )y 

+ (i* - 4 >nw - 4> c + 4>w) £y 


with a corresponding normal gradient, within cell W, 

= (i c - i w ) + (i N - i m ~i c + i w ) v < 24 ) 


Using local, cell-centered coordinates, similar formulae hold for each cell; in particular, cell 
SW formulae can be obtained from (23) and (24) by shifting all indexes “south” by one unit 
in (16), (23), and (24). 

As before, the convective flux integral is split into three geometrically distinct parts. 
In this case. 


When rearranged as 


0.5 r 0.5 -| 

j [ | <t>(^,v) dy\ 


l = 


<HS,V) dr, d£ 


d£ 


= c 


0.5 -c„ -0.5 

0 5 r- 

J ^ + ( <t> c - 

0.5 -c„ 

i w + - K) 


A C xw 


+ K) - -j- (<t>c - K)} 


(25) 


(26) 

(27) 

(28) 


this will be recognised as the one-dimensional Lax-Wendroff flux [3]. The second (negative) 
contribution over the triangular area in cell W is 



0.5 f 0.5 



0-5-c„ ,„({) 


<t>(Z,v) dr, d£ 


(29) 
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Using (17), this becomes, after some rearrangement, 


h = 


c c 

xw yw 


J (<t>c + <t>W + + 4w) 


- ^ (4>N + &C - <t>NW - <t>w) 


■J2 &N + NW ~ 4>c ~ <t>w) 


™ ^ (<t> N - 4 >mw - 4>c + 4> w ) 
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( 30 ) 


Then, I 3 is obtained from I 2 by shifting all indices south by one unit. 
The diffusive flux is computed in a similar way. In particular 




(it) = 

- — f [ 

' di; 



(31) 


FIP 


This is also conveniently split into three parts. In this case 

“ 7^ h = ~ “w dc ~ K) (32) 

C xw 

which will be recognised as the classical, second-order, one-dimensional expression for the 
diffusive flux across the west face. But there are also contributions from transverse 
convective coupling. In particular, 


a 

W J _ 

A = QL 


-J- + 4> C - <t>NW ~ <t>w) 


C 2 - - - - 

- (&N - 4>NW - 4>C + 4>w) 


(33) 


and the corresponding / 3 term is again obtained by shifting all indices south by one unit. 
The total west-face convective-plus-diffusive flux is then 
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BILINEAR DOWNWIND: 


FLUX„(iJ) 


\ tic - C -f ti c - K) 

- - *s) * ti m - *»,)] 

* c -^[ti K - ^ - ti m - «] 

+ IT ” 2 ^ C + ^s) + (4>MV ~ 2i w + 05w)] 

_ [(^ - 24> c + is) - dw - 24> w + isJ] 


(<t>c ~ K) ~ C -f [(** ~ *s) - (<*W - <M] 

+ c -f [d N - 2i c + - (* w - 2i w * iji 


(34) 


The interesting thing about this formula is that (referring to Figure 4) every term is face- 
centered in both x and y directions. In this case, the downwind-weighted sub-cell 
interpolation is “balanced” by the natural upwinding involved in the flux-integral 
calculation. The resulting convective-plus-diffusive flux is independent of velocity direction 
— just like the Lax-Wendroff method in one dimension. The overall update equation 
involves the square, nine-point stencil, centered on C. For pure convection at constant 
velocity, the update is identical to that of a semi-Lagrangian scheme using the same nine- 
point stencil. 

Although the semi-Lagrangian convection scheme can be obtained from the flux- 
integral form, the reverse is not true. This is easily seen by writing out the complete update 
based on (34). Notice how the c x {c/: y ) term from the east- west flux difference combines with 
the c y (c x ) term from the north-south flux difference. The purely convective von Neumann 
stability region is again the square, given by (14). Stability regions in the (c x , c y ) plane for 
finite values of a are discussed elsewhere [4]. 
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UNIFORMLY THIRD-ORDER POLYNOMIAL INTERPOLATION ALGORITHM 


Assuming symmetrically weighted cellwise quadratic sub-cell interpolation within each cell 
leads to a uniformly third-order polynomial interpolation algorithm for convection and 
diffusion (under constant-coefficient conditions). In a variable (but solenoidal) convecting 
velocity field, with possibly variable diffusivity, the algorithm is no longer formally third- 
order accurate; however, the practical accuracy is significantly better than that of formally 
second-order schemes. Phase accuracy, in particular, is excellent — just as in the case of the 
corresponding one-dimensional QUICKEST scheme [5,6,7]. 

Within cell W, the quadratic interpolation takes the (velocity-direction-independent) 

form: 

= 4>w - dc + 4>SW + $ww + “ 4w) 

+ ~ + ~ 2<l> w + 0ww j £2 

+ AW ~ + I * «" ~ + (35) 

Note that this satisfies the cell-average constraint, (8). Also note that the nodal value, <f> w , 
is not the same as the cell average; in fact, 

<t>w = <K0,0) = 4 W - d c + 4>sw + 4>ww + 4>nw - 4w) < 36 ) 


The x-direction gradient within cell W is 


d<f> 


- c - ~ - (J c - 2i w 


+ 



(37) 


Convective Flux in the Absence of Diffusion 

For pure convection, the west-face flux is computed in the usual way. As perhaps, by now, 
expected, the first component of the flux integral generates the one-dimensional (in this case, 
QUICKEST formula: 

2 

A = (♦. + <♦« - K) - (i - (*c - 2K * *J[ (38) 
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Note the appearance of the upwind-weighted normal-curvature term (resulting from the 
natural upwinding inherent in the flux integral). Integration over the respective triangular 
areas in cells W and SW introduces several cross-difference terms. The final form of the 
purely convective flux is (referring to Figure 4): 

QUADRATIC (CONVECTION): 


FLU X w (j,j) 


\ (^c + <t> w ) - (<t> c - 4> w ) 

- (i^)(5c - 2i w + ~4> ww ) 

- - *sw) 


(I - £fn)(4> 
'4 3 ' 

c & w $ s + & sw ) 

'(? - %)<■+** - « 

1 1 

(4> c - 2 4> w + 4 > ww ) 

*12 IP 

- ~ (^5 “ ^SW + ^SWw)- 


2 

+ Cyw ^12 ~ ’24^^ NW ~ ^ w + ~ tssw) 


(39) 


It is instructive to identify the role played by each of these terms. The first three 
terms represent the one-dimensional (QUICKEST) contribution; note that all the remaining 
terms contain a coefficient. The fourth term in (39) is the transverse gradient that 
previously appeared in the first-order scheme. This is, in fact, a second-order term — just 
as is the normal gradient (second term). The fifth term represents twist, an interaction 
between normal and transverse convection. The next term is a transverse-curvature 
contribution. The final two terms in (39) are actually fourth-order contributions. Dropping 
these terms does not affect the formal (constant-coefficient) third-order accuracy of the 
overall update equation. However, they do affect the stability of the scheme. Without the 
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higher-order terms, the purely convective stability region is approximately the diamond- 
shaped region 


Including these terms (that arise naturally in the flux-integral formulation) results once again 
in the square stability region given by (14). 

Third-Order Diffusive Flux 

Applying the usual three-part-integral procedure to (37) results in the following diffusive 
flux: 


The first term is the classical (one-dimensional) second-order first-difference across the face. 
The second, normal-curvature term, represents the effect of normal convection on the time- 
averaged normal gradient; this is a third-order convection-diffusion cross-coupling term that 
also appears in the one-dimensional QUICKEST formula [5]. The third, twist, term 
represents the coupling effect of transverse convection; this is a (two-dimensional) third- 
order term. The final term is actually a (partial) fourth-order cross-coupling term, kept, 
once again, because of enhanced stability properties [4]. 

Diffusive Contribution to the Convective Flux 

The convection-diffusion coupling terms just described represent the effects of convection 
in estimating the diffusive flux. They arise naturally in the flux-integral formulation. For 
uniformly third-order consistency, one also needs to estimate the analogous cross-coupling 
effects of diffusion on the convective flux. Because of the assumed curvature in the sub-cell 


Kl + \c y \ < 1 


(40) 


(4>C - <t>w) - ~Y (4>c ~ 2 4> w + <?w) 
~ ).vg = ” ~ ~Y iic " ~ is + 4>sw) 



(<t> c 20 ^, + 


(41) 
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interpolation, diffusion changes the value of 4> over At as it is being convected through a 
particular cell face. As shown in [8], this change is given by Performing the 

usual three-part flux-integral calculation leads to an additional diffusive contribution to the 
convective flux of the form 


~Y [(^c “ 2< t>w + ^W) + (<t>Nw ~ 2 $w + tsw)] 


+ 


a c 

w yw 


(<f>C 2</> w + <t > ww ) ( 4>s + $ sww ) 

- + (^W ~~ + ^sw ~ & ssw ) 


(42) 


which must be added to (39) to give the total convective flux. The convective-plus-diffusive 
flux at the west face is thus the sum of (39), (41), and (42). A von Neumann stability 
analysis of the constant-coefficient overall update algorithm [4] shows that the useful region 
in ( c x , c y , a) space is given, as a minimum, by the “cylinder”: 

|c| < 1, 0 < a < 0.25 ( 43 ) 


Reference [4] also describes a simple and inexpensive (vectorisable) upwinding strategy 
based on the “generic stencil” technique, using FORTRAN functions NINT and SIGN. 


CONVECTIVE TEST PROBLEM 

The three schemes described here have been applied to the rotating Gaussian hill problem 
under purely convective conditions (a=0). Details of the mesh, time-step, and other 
parameters are given in Reference [4], Figure 6 shows the initial state in a superfine-grid 
rendering. The cell-average initial data is shown in Figure 7. Figure 8 shows the first-order 
results after one (anticlockwise) rotation. The “Lax-Wendroff” cell-average results are 
shown in Figure 9, with a close-up of the peak region, showing the sub-cell behaviour, in 
Figure 10. Figures 11 and 12 give the corresponding UTOPIA results. As mentioned 
before, the first-order scheme is far too artificially diffusive to even be considered for 
practical application. As in one dimension, the Lax-Wendroff-type scheme is excessively 
dispersive, showing significant phase-lag errors in the “wake”. By contrast, UTOPIA has 
good accuracy and excellent phase behaviour, just as in the one-dimensional case [5,7]. 
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Extensive other studies of purely convective and convective-diffusive test-problems 
at a number of grid-refinements have consistently shown the superiority of UTOPIA. Not 
surprisingly, it is more accurate than lower-order schemes. It is also more expensive — per 
mesh point calculation. The important conclusion, however, is that, for a prescribed 
accuracy , UTOPIA can be used on a much coarser mesh (with a concomitantly larger time- 
step) — the overall cost is then much lower than that of lower-order schemes. 

CONCLUSION 

The flux-integral method is a powerful technique for estimating convective-plus-diffusive 
fluxes in a strictly conservative formulation of an explicit, single-step update formula for the 
multidimensional convection-diffusion equation. The assumption of locally constant 
convecting velocities near each control-volume face means that the formal accuracy of the 
convection terms is at most second order — e.g., the convecting velocity field could be 
staggered in time by Ar/2 with respect to transported scalars. Variable diffusivity would 
typically be lagged, and therefore only first-order accurate in time. However, spatial 
accuracy is enhanced by using higher-order sub-cell interpolation. Under constant- 
coefficient conditions, an Nth-order sub-cell interpolation leads to an ( N+ l)th-order accurate 
convection-diffusion scheme — in both space and time. 

Cellwise constant sub-cell interpolation leads to an unworkable (artificially diffusive) 
first-order convection scheme. Cellwise linear or bilinear interpolants generate second-order 
convection-diffusion schemes. Downwind-weighted bilinear interpolation gives a multi- 
dimensional analogue of the Lax-Wendroff scheme. However, because of the unsymmetrical 
weighting of the interpolant, this leads to a highly dispersive convection scheme with strong 
phase-lag errors — just as in one dimension. Cellwise quadratic interpolation (independent 
of velocity direction) leads to a very accurate convection-diffusion scheme with excellent 
phase behaviour. Under constant coefficient conditions, this is a uniformly third-order 
polynomial interpolation algorithm (UTOPIA). Since highly accurate solutions can be 
obtained on relatively coarse grids, UTOPIA is much more cost-effective than lower-order 
methods. 
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Because fluxes are estimated directly, (for a given sub-cell interpolation) the flux- 
integral method produces unique formulae for the fluxes. This is in contrast to other 
methods that have been used to construct conservative multidimensional convection (or 
convection-diffusion) schemes. For example, Ekebjaerg and Justesen [9] developed a 
nominally third~-order, two-dimensional convection-diffusion scheme by successive 
elimination of truncation error terms arising from a lower-order scheme. The 
nonconservative single-step explicit update was then rewritten in a conservative pseudoflux- 
difference form . But the pseudofluxes chosen by Ekebjaerg and Justesen are not unique; 
except for first-order (and some simple second-order) schemes, there is, in general, no 
unique way of rewriting a nonconservative update in a conservative flux-difference form [4], 

Recently, Rasch [10] has used, as a starting point, a (constant-coefficient) third-order 
semi-Lagrangian convection scheme, aind then rewritten the update in conservative form. 
Recognising the nonuniqueness problem, Rasch uses weighting parameters to generate a 
family of possible pseudoflux-difference algorithms. For certain choices of the weights, the 
(purely convective form of the) Ekebjaerg-and -Justesen scheme can be retrieved. For other 
weights, Rasch’s convection scheme is equivalent to (the convective part of) that used by the 
present authors in Reference 8; this is also the form used by Rasch. For still other choices 
of the weights, a nominally third-order convective flux equivalent to (39) — but with the last 
two (higher-order) terms removed — can be obtained. The conservatively rewritten semi- 
Lagrangian approach does not generate diffusive fluxes. 

In a very recent manuscript, LeVeque [11] has used a technique for purely convective 
flows similar, in some respects, to the flux-integral method described here. For a constant 
convecting velocity field, LeVeque’s Method V is a ten-point third-order two-dimensional 
convection scheme. This is the minimum number of points needed for third-order accuracy. 
The overall convective update is equivalent to the semi-Lagrangian scheme used as the 
starting point for Rasch’s method; it is also equivalent to that of Ekebjaerg and Justesen. 
LeVeque’s Method VI is equivalent to the purely convective portion of the flux-integral 
method developed here; i.e., (39). LeVeque does not consider diffusive fluxes. 
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For simplicity, the present paper has been confined to two dimensions. It should be 
clear that the flux-integral method generalises to three dimensions in a straight-forward 
manner. The technique can be used with even higher-order sub-cell interpolation 
(presumably using velocity-direction-independent interpolants in order to minimise 
dispersion); however, higher-order diffusive-diffusive and convective-diffusive cross- 
coupling terms appear to be quite complex. Further research is needed to extend the method 
to larger time steps. Finally, it should be pointed out that shape preservation in the sub-cell 
interpolation automatically results in a positivity-preserving conservative multidimensional 
formulation — thus obviating the need for ad hoc flux-limiter constraints. This is an area of 
current research. 
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Figure 1 


Schematic of particle paths flowing into the west face of cell C. 



Figure 2 The flux-integral parallelogram in the vicinity of the west face of cell C; 

Cxw > Cyw ^ 0 * 
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Figure 3 Flux-integral parallelograms on each face of cell C. 
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Figure 4 Twelve cells in the vicinity of the west face of cell C; c x , c y > 0. 
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Figure 5 



Definition of local normalised coordinates within cell W. 
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Figure 6 Initial Gaussian distribution shown on a very fine grid. This is also the exact 
(purely convective) solution after an integral number of rotations; 

= 1 - 000 . 
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Figure 7 
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Initial state of cell-average values, shown as a two-dimensional histogram on 
the computational grid; ^ = 0.991. 
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Figure 8 Solution after one (anticlockwise) rotation for the first-order method. In this 
case, the histogram of cell-average values also represents the cellwise constant 
sub-cell interpolation. = 0.152; = 0.0. 
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Figure 9 Solution of cell-average values after one rotation using the second-order 
convection scheme; 4 >„„ = 0.787; 4^ = -0.149. 
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Figure 10 Close-up of the downwind-weighted bilinear sub-cell interpolation in the peak 
region, showing discontinuities across cell faces; 4>„. r = 0.823; 

^ = -0.176. 
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Figure 11 


Solution of cell-average values after one rotation using the purely convective 
form of UTOPIA; = 0.804; = -0.008. 
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Figure 12 Close-up of the quadratic sub-cell interpolation in the peak region, showing 
(small) discontinuities across cell faces; 4 ^ = 0.810; = -0.010 (not 

shown). 
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